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1.  INTRODUCTION 


Increasingly,  today1 s  manufacturing  companies — both  large 
and  small — are  finding  it  necessary  to  adopt  computer-based  planning 
systems  to  monitor  their  production,  distribution,  and  inventory 


XPDI)  functions 


jK 


This  paper  describes  an  optimization  model 


that  captures  the  multi-period  aspect  of  many  PDI  planning  systems. 
This  model  is  based  on  the  well-known  pure  network  flow  (transship¬ 


ment)  problem 


The  model  can  be  used  to  optimally  select : 


^ro^ctjLon^levels  at  multiplexing  ^ 

(primary)  distribution  patterns  between  plants  and  storage 


facilities,^ 

J&i  inventory  ^Levels  at  multiple  storage  sites , 

^(secondary)  distribution  pattems  between  storage  facilities^ 
»  and  customer  demand  points, 

^customer  satisfaction  levels.  ^ 

^  /1\ 

The  model  determines  the  optimal  solution  for  each  of  these  factors 
in  each  time  period.  By  simultaneously  considering  all  time  periods 
(e.g.,  months)  within  the  span  of  the  planning  horizon  (e.g,,  fiscal 
year)  the  model  is  able  to  globally  optimize  the  solution.  That  is, 
the  model  can  tailor  production  and  inventory  schedules  to  account 


for  fluctuating  demands  and  prices  over  the  planning  horizon! 
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Since  this  PDI  planning  model  is  based  on  the  pure  network 
flow  problem ,  it  can  be  solved  using  any  of  the  numerous  capacitated 
transshipment  computer  codes  such  as  PNET  [  1,  10,  11],  GNET  [6  ], 
RNET  [13],  or  SUPERK  [3  ].  However,  in  terms  of  computer  memory  re¬ 
quirements,  each  of  these  codes  suffers  considerably  when  the  number 
of  time  periods  is  moderately  large.  This  is  due  to  the  fact  that 
multi-period  planning  models  generally  require  a  separate  "copy*1  of 
the  basic  PDI  network  for  each  time  period. 

A  specialized  implementation  of  the  primal  simplex  network 
code  PNET  has  been  developed  that  overcomes  this  computer  memory 
handicap.  This  implementation,  PNET-MP,  is  unique  in  that  it  only 
stores  a  single  copy  of  the  basic  PDI  network,  thus  minimizing  the 
impact  of  a  planning  horizon  with  multiple  time  periods.  PNET-MP 
generates  the  complete  multi-period  network  as  it  is  needed.  Many 
of  the  efficient  list  structures  that  have  been  incorporated  in  the 
other  primal  simplex  codes,  such  as  PNET,  GNET,  and  RNET,  were 
adopted  in  this  multi-period  code. 

PNET-MP  was  developed  solely  to  reduce  the  computer  memory 
requirements  for  this  important  class  of  planning  problems.  A  de¬ 
gradation  in  solution  speed,  relative  to  PNET,  was  anticipated  since 
extra  data  manipulations  must  be  carried  out  in  order  to  generate 
the  multiple  copies  of  the  PDI  network  from  the  single  copy  that  is 
actually  stored.  However,  as  will  be  illustrated  in  Section  6, 
PNET-MP  is  actually  more  efficient  than  its  general-purpose  parent 
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code  for  this  special  class  of  planning  problems*  This  is  due  to 
the  modifications  of  the  pivot  strategy  that  were  mandated  by  the 
code’s  restricted  access  to  the  complete  multi-period  network. 

2.  PROBLEM  DESCRIPTION 

In  this  section  a  general  statement  of  the  specific 
problem  class  under  consideration  is  presented.  This  type  of  problem 
arises  in  many  production,  distribution,  and  inventory  planning 
systems.  In  general,  the  problem  is  to  determine  a  minimum  cost 
routing  of  units  of  a  homogenous  commodity  through  a  distribution 
network  over  a  (finite)  multi-period  planning  horizon.  The  specific 
units  being  routed  may  be  oil  tankers,  freight  cars,  aircraft,  or 
military  personnel.  It  must  be  assumed  that  the  units  are  suffi¬ 
ciently  similar  to  be  treated  within  a  single  commodity  framework. 

For  instance,  the  units  may  actually  be  forty  foot  refrigerated  box 
cars  or  second  lieutenants  with  infantry  skills. 

Any  distribution  network  consists  of  a  set  of  nodes  N  and 
a  set  of  arcs  A  (see  Figure  1).  Each  node  i  £  N  of  a  distribution 
network  represents  a  specific  geographic  location  such  as  a  production 
facility,  a  warehouse,  a  customer,  or  a  military  location.  Each  arc 
k  e  A  represents  an  allowable  distribution  link  between  an  ordered 
pair  of  nodes  That  arc  k  represents  a  link  from  node 

i^  e  N  to  node  j  e  N.  Arcs  can  correspond  to  shipping  channels,  rail* 
road  tracks,  flight  path  corridors,  or  personnel  reassignments. 

Each  node  of  the  distribution  network  has  a  schedule  of 


FIGURE  1 

SIMPLE  DISTRIBUTION  NETWORK 


supplies  and  demands  for  units  at  specific  points  In  time  over  the 
planning  horizon.  Some  or  all  of  the  nodes  in  the  network  can  be 
used  to  store  units  until  they  are  needed.  These  storage  nodes  have 
both  a  marginal  storage  cost  and  a  capacity  limitation  associated 
with  them.  For  instance,  in  Figure  1,  the  node  representing  the 
Dallas  warehouse  could  have  a  storage  cost  of  $10  per  unit  per  period 
and  a  storage  capacity  of  5  units  per  period.  Each  arc  of  the  net¬ 
work  has  a  marginal  distribution  cost  as  well  as  a  required  tra- 


versal  time.  All  arcs  are  assumed  to  be  uncapacitated.  In  Figure  1, 
the  arc  from  the  Dallas  warehouse  node  to  the  Houston  customer  node 
might  have  a  marginal  cost  of  $50  per  unit  and  a  traversal  time  of 
ten  hours. 

The  only  fundamental  difference  between  this  problem  and 
the  standard  pure  network  problem  is  that  the  arcs  have  a  delay  or 
traversal  time  associated  with  them.  This  places  a  restriction  on 
the  movement  of  units  over  time.  This  restriction  is  not  captured 
by  the  standard  pure  network  model. 

In  order  to  simplify  the  problem  analysis,  the  planning 
horizon  will  be  partitioned  into  T  equal  length  time  periods.  De¬ 
pending  upon  the  characteristics  of  the  specific  problem  being 
solved,  these  periods  may  be  hourly,  daily,  or  even  weekly. 

The  supplies  and  demands  for  units  at  each  node  in  the 
network  are  assumed  to  take  place  at  the  end  of  the  time  periods. 

For  instance,  two  empty  freight  cars  could  become  available  Tuesday 
morning  in  Chicago  (supply),  and  one  empty  freight  car  could  be  re¬ 
quired  by  Friday  morning  in  Houston  (demand).  From  period  to  period, 
the  total  number  of  available  units  may  fluctuate  as  they  are  being 
added  to  and  removed  from  the  distribution  system. 

In  order  to  simplify  the  problem  analysis,  it  will  be  as¬ 
sumed  that  all  units  must  arrive  at  a  node  by  the  end  of  each  period. 
That  is,  a  unit  cannot  be  left  stranded  between  nodes  at  the  end  of  a 
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period.  In  many  cases  this  Is  a  realistic  assumption.  In  those 
cases  where  this  assumption  Is  not  reasonable ,  long  arcs  can  be  broken 
up  by  introducing  artificial  nodes. 

Many  real-world  distribution  systems  cover  large  geo¬ 
graphic  regions.  This  places  severe  limitations  on  the  ability  to 
meet  the  demand  for  units.  For  example,  it  may  not  be  possible  to 
use  the  empty  freight  cars  in  Chicago  to  meet  the  demand  in  Houston 
even  though  they  are  available  three  days  in  advance.  This  model 
seeks  to  capture  this  restriction  on  the  length  of  allowable  shipping 
routes  through  a  distribution  system. 

Two  different  formulations  of  this  problem  are  presented  in 
the  following  sections.  The  first  is  an  integer  linear  programming 
formulation  of  a  multicommodity  network  model.  The  second  is  a 
capacitated  transshipment  formulation  of  a  single  commodity  network 
model.  Technically,  both  formulations  are  equivalent.  However,  in 
terms  of  computational  feasibility,  the  second  formulation  is  over- 
whemingly  superior. 

3.  MULTICOMMODITY  FORMULATION 

An  integer  linear  programming  formulation  of  the  problem  is 
presented  in  this  section.  The  formulation  is  essentially  an  integer 
multicommodity  model.  The  complexity  of  this  approach  is  a  direct 
result  of  the  stipulation  that  the  total  traversal  time  of  each  ship¬ 
ping  route  cannot  exceed  the  length  of  the  time  period. 

This  approach  is  a  multicommodity  model  because  it  is 
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necessary  to  keep  track  of  the  origin  and  destination  of  each  unit 
as  it  is  moved  through  the  network.  During  any  given  period,  the 
units  that  are  being  shipped  across  an  arc  may  have  originated  (at 
the  beginning  of  the  period)  at  any  of  the  nodes  in  the  network, 
and  may  be  destined  (at  the  end  of  the  period)  for  any  of  the  other 
nodes.  However,  due  to  the  traversal  time  restriction,  many  com¬ 
binations  of  origin  and  destination  nodes  are  not  feasible.  For 
instance,  it  may  be  physically  impossible  to  move  an  empty  freight 
car  from  Chicago  to  Houston  in  a  single  day. 

In  order  to  differentiate  between  the  various  units  tra¬ 
versing  an  arc,  an  integer  flow  variable  must  be  introduced. 

The  index  k  denotes  the  arc  number,  i  is  the  original  origin  node, 
j  is  the  final  destination  node,  and  t  is  the  time  period.  That  is, 
i  is  the  node  where  the  unit  was  located  at  the  start  of  period  t, 
and  j  is  the  node  where  the  unit  will  be  located  at  the  end  of  period 
t.  It  should  be  noted  that  nodes  i  and  j  are  not  necessarily  the 
endpoints  (ifc  and  j  )  of  arc  k.  (See  Figure  2.)  That  is,  a  unit 
may  be  shipped  from  node  i  to  node  j  via  many  arcs  during  the  same 


FIGURE  2 
SHIPPING  ROUTE 
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period  provided  that  the  total  traversal  time  of  the  shipping  route 
does  not  exceed  the  length  of  the  period. 

The  storage  capability  of  a  node  can  be  modeled  as  a 
loop  arc  from  the  node  to  itself.  Then  the  number  of  units  stored 
at  node  v  during  period  t  is  given  by  where  k  is  the  index 

of  the  loop  arc. 

Each  arc  has  a  marginal  cost  associated  with  it.  For 
the  distribution  arcs  (i^  f  j^)  this  is  the  marginal  shipping  cost, 
and  for  the  loop  arcs  (i^  =  j^)  this  is  the  marginal  storage  cost. 
For  some  problems,  these  marginal  costs  may  vary  over  the  planning 
horizon.  In  these  cases  it  is  necessary  to  introduce  a  time  index, 


Each  loop  arc  has  an  upper  bound  corresponding  to  the 
storage  capacity  of  the  associated  node.  The  upper  bound  for 
is  given  by  or  if  the  storage  capacity  varies  over  time. 

Each  distribution  arc  has  a  traversal  time  L  associated 

k 

with  it.  It  is  convenient  to  scale  these  times  relative  to  the 
length  of  a  period.  For  instance,  if  a  daily  period  is  used  and 
arc  k  requires  six  hours  to  traverse,  then  “  .25. 

Each  node  has  a  supply  or  demand  for  units  at  the  end  of 
each  period.  Let  b^  be  node  v’s  net  demand  at  the  end  of  period  t. 

It  is  notationally  convenient  to  introduce  two  arc  index 
sets  for  each  node.  The  first  set 


A(v)  -  {k  |  ik  -  v,  Jk  +  v} 
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is  the  set  of  distribution  arcs  that  originate  at  node  v,  and  the 


second  set 


B(v)  =  {k  |  ifc  t  v,  Jk  =  v} 

is  the  set  of  distribution  arcs  that  terminate  at  node  v.  These 
sets  do  not  include  the  loop  arcs. 

In  addition  to  the  integer  flow  variable  ,  a  binary 

variable  is  needed  in  order  to  capture  the  traversal  time 

restriction  on  feasible  shipping  routes.  This  variable  is  used 
to  identify  the  arcs  that  are  used  to  ship  units  from  node  i  to 
node  j  during  period  t. 

The  T  time  periods  are  numbered  from  1  to  T,  and  the 
T  +  1  points  in  time  at  the  ends  of  periods  are  numbered  from  0 
to  T.  That  is,  the  point  in  time  at  the  start  of  period  t  is 
numbered  t  -  1,  and  the  point  in  time  at  the  end  of  period  t  is 
numbered  t. 

The  set  of  network  nodes  is  denoted  by  N.  The  set  of 
distribution  arcs  is  denoted  by  A. 

The  objective  of  this  problem  is  to  minimize  the  total 
cost  of  distribution  and  storage  over  the  T  time  periods.  This 
objective  can  be  stated  as 


Minimize  J  J  ck  1  1  \±i  +  I  2  ckxkii 

t=l  k£A  K  ieN  jeN-i  J  t=l  icW  K 


where  the  first  term  is  the  total  distribution  cost  and  the  second 
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term  is  the  total  storage  cost.  The  quantity 

ieN  jeN-i 


is  the  total  number  of  units  shipped  across  arc  k  during  period  t. 

There  are  five  basic  types  of  constraints  for  this  integer 
linear  programming  formulation  of  the  problem.  These  constraint 
types  are: 

!•  Conservation  of  flow  originating  or  terminating  at  node  v  at 
time  t, 

II.  Conservation  of  flow  passing  through  node  v  during  period  t, 

III.  Flow  variable/binary  variable  relationships, 

IV.  Shipping  route  total  traversal  time  limitations,  and 

V.  Node  storage  capacities. 

Each  constraint  type  will  be  presented  separately. 


Type  I  Constraints 

These  conservation  of  flow  constraints  affect  only  the 
units  that  start  a  period  or  end  a  period  at  node  v.  These  con¬ 
straints  guarantee  that  the  number  of  units  that  terminate  at 
node  v  at  the  end  of  period  t  equals  the  number  of  units  that 
originate  at  node  v  at  the  start  of  period  t  +  1  plus  the  net  de¬ 
mand  for  units  at  node  v  at  the  end  of  period  t.  There  is  one 
type  I  constraint  for  each  combination  of  network  node  and  point 
in  time.  The  general  form  of  the  constraint  is 
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Type  III  Constraints 


These  constraints  enforce  the  proper  relationship  between 
the  flow  variables  an<*  t*ie  binary  variables  Specifically, 

l£  must  be  equal  to  one  unless  is  equal  to  zero.  There  is 

one  Type  III  constraint  for  each  combination  of  distribution  arc, 
origin  node,  destination  node,  and  time  period.  The  constraint 
takes  the  form 

t  <  t 

*kij  *  “Scij 

where  M  is  an  arbitrarily  large  number. 


Type  IV  Constraints 

These  constraints  restrict  the  total  traversal  time  of 
each  shipping  route.  The  validity  of  these  constraints  relies  on 
the  following  observation.  Namely,  since  all  distribution  arcs 
are  uncapacitated,  all  units  being  shipped  from  node  i  to  node  j 
during  period  t  can  travel  along  a  single  (least  cost)  route.  This 
means  that  the  total  traversal  time  of  the  shipping  route  between 
nodes  i  and  j  during  period  t  can  be  measured  by 


1  L 


keA 


k  Stij 


Therefore,  the  restriction  that  the  total  traversal  time  of  the 
shipping  route  between  nodes  i  and  j  is  less  than  the  length 
of  the  period  can  be  enforced  with  the  following  constraint 


.?» **  ■ 1 


keA 


There  is  one  Type  IV  constraint  for  each  combination  of  origin 
node,  destination  node,  and  time  period. 


Type  V  Constraints 

The  storage  capacity  at  each  node  can  be  handled  as  a 
simple  upper  bound  on  the  flow  on  each  loop  arc.  This  constraint 
has  the  form 

x5  =  S 
Tew  v 

There  is  one  Type  V  constraint  for  each  combination  of  network  node 
and  time  period. 

This  formulation,  although  mathematically  correct,  is 
ovbiously  computationally  infeasible.  A  simple  five  period  problem 
with  an  underlying  100  node  and  500  arc  distribution  network  pro¬ 
duces  an  integer  linear  programming  problem  with  nearly  25  million 
constraints  and  50  million  variables.  As  stated  earlier,  the 
complexity  of  this  approach  is  a  direct  result  of  the  restriction 
on  the  total  traversal  time  of  each  shipping  route. 

4.  CAPACITATED  TRANSSHIPMENT  FORMULATION 

In  this  section,  the  problem  introduced  in  Section  2  is 
formulated  as  a  multi-period  capacitated  transshipment  problem. 

This  single  conmodity  formulation  removes  the  majority  of  the  com¬ 
plexity  of  the  multi- commodity  formulation  of  Section  3.  Notably, 
the  integrality  requirements  become  a  natural  consequence  of  the 
unimodularity  property  of  the  network  constraints  [15]. 


This  transshipment  formulation  of  the  probelm  revolves 
around  a  simple  transformation  of  the  underlying  distribution 
network.  The  transformed  network  contains  the  same  nodes  as  the 
original  network,  but  contains  a  different  set  of  arcs.  Roughly 
speaking,  the  transformed  network  has  an  arc  corresponding  to 
each  feasible  shipping  route  in  the  original  distribution  network. 

The  arcs  of  the  transformed  network  can  be  generated 
by  solving  a  set  of  simple  shortest  path  problems.  The  first 


shortest  path  problem  uses  the  traversal  time  as  the  length  of 


each  arc  in  the  original  network.  The  solution  to  this  shortest 
path  problem  provides  the  minimum  traversal  time  of  the  shipping 


route  between  each  pair  of  nodes.  Let  T^  be  the  minimum  traversal 


time  from  node  i  to  node  j .  The  transformed  network  has  a  directed 


arc  from  node  i  to  node  j  if  and  only  if  T^.  5  1.  That  is,  the  trans 


formed  network  has  an  arc  between  each  pair  of  nodes  that  are  less 
than  one  period  apart  in  terms  of  traversal  time. 

The  marginal  cost  of  an  arc  in  the  transformed  network 
is  equal  to  the  sum  of  the  marginal  costs  of  the  arcs  in  the  cor¬ 
responding  shipping  route  in  the  original  network.  If  there  are 
multiple  shipping  routes  between  a  pair  of  nodes  in  the  original 


network,  then  a  shortest  path  problem  with  c^  as  the  length  of 


arc  k  can  be  used  to  determine  the  marginal  cost  of  the  best  ship¬ 
ping  route.  Caution  must  be  exercised  to  make  sure  that  the 
minimum  cost  shipping  route  has  a  total  traversal  time  no  longer 
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than  a  period.  For  instance,  if  three  shipping  routes  between 
nodes  1  and  2  in  the  original  network  have  total  traversal  times 
of  .4,  .9,  and  1.6,  and  have  total  marginal  costs  of  20,  30,  and 
15,  then  the  transformed  network  would  have  a  directed  arc  from 
node  1  to  node  2  with  a  marginal  cost  of  20.  Note  that  the  best 
marginal  cost  (15)  is  associated  with  a  shipping  route  that  re¬ 
quires  more  than  one  time  period  (1.6).  A  transformation  of  this 
type  was  proposed  by  Wagner  [19]. 

Since  all  of  the  distribution  arcs  in  the  original  net¬ 
work  are  uncapacitated,  all  shipping  routes  in  the  transformed  net¬ 
work  are  also  uncapacitated.  Since  the  arcs  in  the  transformed 
network  actually  correspond  to  feasible  shipping  routes  in  the 
original  network,  all  arcs  in  the  transformed  network  are  uncapa¬ 
citated. 

The  transformation  is  completed  by  replacing  the  restric¬ 
tion  that  the  total  traversal  time  of  each  shipping  route  is  less 
than  one  period,  by  the  restriction  that  a  unit  of  flow  can  tra¬ 
verse  at  most  one  arc  per  period  in  the  transformed  problem. 

As  before,  the  storage  capability  at  node  v  is  handled 
by  adding  a  loop  arc  from  node  v  to  itself.  The  marginal  cost  of 
this  loop  arc  is  equal  to  the  marginal  storage  cost  at  node  v. 

This  arc  has  an  upper  bound  equal  to  the  storage  capacity  of  node 
v.  This  expanded  transformed  network  is  called  the  basic  network . 

Let  N  be  the  set  of  nodes,  and  A  be  the  set  of  arcs  in 
the  basic  network. 
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The  number  of  units  shipped  across  arc  k  during  period 
t  is  denoted  by  x^.  The  marginal  cost  of  arc  k  is  given  by  c^ 


(or  c  ) ,  and  the  upper  bound  on  the  flow  on  arc  k  is  given  by  u 


(or  .  So  if  arc  k  corresponds  to  the  loop  arc  for 

node  v.  Otherwise,  ®  <». 

The  net  demand  for  units  at  node  v  at  the  end  of  period 

t  is  denoted  by  bt . 

v 

As  before,  it  is  convenient  to  introduce  two  arc  index 
sets.  The  first  set 


A(v)  -  {k  |  ik  =  v} 

is  the  set  of  arcs  that  originate  at  node  v  in  the  basic  network, 
and  the  second  set 

B(v)  -  {k  |  jfc  -  v} 


is  the  set  of  arcs  that  terminate  at  node  v  in  the  basic  network. 
The  loop  arcs  are  included  in  these  sets. 

The  objective  function  for  this  formulation  of  the 
problem  is  given  by 


Minimize 


T 

1 

k«l 


k£A 


t 

Vk 


This  formulation  has  two  basic  constraint  types : 

A.  Conservation  of  flow  at  node  v  at  time  t,  and 

B.  Arc  capacities. 

These  constraints  are  considered  Individually. 
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Since  this  formulation  requires  that  the  traversal  of 
any  arc  of  the  basic  network  consumes  a  complete  time  period,  no 
units  will  pass  through  an  intermediate  node  during  a  period. 
This  eliminates  the  need  for  a  constraint  analogous  to  the  Type 
II  constraint  of  Section  3.  All  conservation  of  flow  is  handled 
by  the  following  version  of  the  Type  I  constraint  of  Section  3. 
There  is  one  Type  A  constraint  for  each  combination  of  network 
node  and  point  in  time.  The  general  form  of  the  constraint  is 


The  storage  capability  at  the  nodes  can  be  handled  as  a 
simple  upper  bound  constraint  on  the  loop  arcs.  There  is  one  Type 
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B  constraint  for  each  combination  of  network  node  and  time  period. 
The  constraint  has  the  form 


xk  -  Uk 


The  complete  algebraic  statement  of  the  capacitated  trans¬ 
shipment  formulation  of  the  problem  is  given  as 


Minimize  £  £ 


t~l  keA 


:k\ 


subject  to: 


-I 

k£A(v) 


1 

*k 


I  *5-1 


k£B(v) 


k£A(v) 


=  b  for  all  v  £  N 
v 


*t+1  =  bC  for  all  v  £  N  and 
k  v 


1,2 . T  -  1 


ksB(v) 

°**k*"k 


b  for  all  v  £  N 
v 


for  k  £  A  and 
t  *  1,  2,  ....  T 


This  problem  is  easily  seen  to  be  a  capacitated  pure  transshipment 
problem  if  it  is  restated  in  matrix  notation.  In  order  to  do  this, 
two  m  x  n  matrices  (F  and  T)  will  be  introduced.  Essentially, 
they  represent  the  origin  and  destination  nodes  of  the  arcs  of  the 
m  node  and  n  arc  basic  network.  The  element  in  the  vC^  row  and  kC^ 
column  of  the  origin  matrix  F  is  defined  by 

-1  if  i,  -  v 
k 

0  otherwise 
th  th 

and  the  element  in  the  v  row  and  k  column  of  the  destination 
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matrix  T  is  defined  by 

1  if  Jk  -  V 

0  otherwise 

So  each  column  of  F  and  T  has  a  single  non-zero  element. 

t  t  t  t  T 

Let  X  be  the  vector  (x^,  x^,  ...»  x^)  of  arc  flow 

during  period  t.  Let  F^  and  Ty  denote  the  vth  row  of  F  and  T, 

respectively.  Then 

F  xc  =  -  7  x*; 

v  keA(v) 

and 

T  Xt  -  2  x} 

V  keB(v)  k 

Let  C  be  the  vector  (c^  c2 . cn),  U  be  the  vector  (u^  u2,  ..., 

X  t  t  t  t  T 

u  )  ,  and  b  be  the  vector  (b. ,  b0 ,  ...,  b  )  .  Then  the  capacitated 
n  14  m 

transshipment  problem  can  be  restated  in  matrix  terms  as 

Minimize  CX1  +  C X2  +  CX3  . . .  +  CXT_1  +  CXT 
subject  to:  ^1  m  b0 

TX1  +  FX2  -  b1 

2  3  2 

TX  +  FX  -  b 


T-l  T  T-l 
TX  +  FX  «  b 

T  T 
TX  -  b 


0  2  XC  £  U  for  t  -  1,  2 . T 


.■J 
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Ignoring  the  simple  upper  bound  restrictions  on  the  flow 
variables,  each  column  of  the  constraint  matrix  of  this  problem 
has  two  non-zero  elements.  Column  k  +  (t  -  l)n  is  associated 
with  variable  This  column  has  a  -1  element  in  row  i^  + 

(t  -  l)m,  and  a  4*1  element  in  row  +  tm.  Therefore,  this  prob¬ 
lem  can  be  seen  to  be  a  capacitated  pure  transshipment  problem 
with  M  *  (T  +  l)m  nodes  and  N  *  Tn  arcs.  This  problem  will  be  re¬ 
ferred  to  as  the  master  network  problem .  The  nodes  and  arcs  of 
the  master  network  problem  are  called  master  nodes  and  master 
arcs .  The  master  network  essentially  consists  of  T  copies  of  the 
original  m  node  and  n  arc  basic  network.  One  master  network  node 
is  generated  for  each  combination  of  basic  network  node  and  point 
in  time,  and  one  master  network  arc  is  generated  for  each  combina¬ 
tion  of  basic  network  arc  and  time  period. 

Clearly  this  formulation  of  the  original  problem  is 
superior  to  the  integer  multi-commodity  formulation  in  Section  3. 

A  direct  comparison  of  the  relative  sizes  of  the  two  problem  formu¬ 
lations  is  not  particularly  easy,  since  the  capacitated  transship¬ 
ment  formulation  Is  based  on  a  transformed  network.  The  number  of 


arcs  in  the  transformed  network,  or  basic  network,  is  a  function  of 


the  number  of  feasible  shipping  routes  in  the  original  distribution 
network.  To  simplify  the  analysis  of  the  relative  problem  sizes, 
let  0  denote  the  average  number  of  arcs  that  originate  at  a  node 
in  the  basic  network.  That  is,  0  is  the  average  number  of  nodes 
that  can  be  reached  within  one  time  period  from  a  given  node. 


In  Section  3,  it  was  shown  that  a  five  period  problem 
with  a  100  node  and  500  arc  distribution  system  generated  an 
integer  linear  programming  problem  with  nearly  25  million  constraints 
and  50  million  variables.  If  0  *  10,  then  this  same  problem  can  be 
formulated  as  a  capacitated  pure  transshipment  problem  with  a  600 
node  and  5,000  arc  master  network.  Even  if  0  is  as  large  as  fifty, 
the  master  network  only  has  600  nodes  and  25,000  arcs.  Clearly,  the 
single  comnodity  formulation  of  the  problem  can  be  readily  solved 
using  any  of  the  available  capacitated  transshipment  computer  codes. 

A  small  example  is  given  in  the  next  section.  In  Section 
6,  a  specialization  of  the  efficient  primal  simplex  network  code 
PNET  is  presented.  This  specialized  code  PNET-MP  can  be  used  to 
solve  much  larger  multi-period  problems  than  can  be  solved  by  PNET 
or  any  of  the  other  in-core  network  codes.  This  specialized  code 
is  extremely  valuable  when  either  T  or  0  is  large. 

5.  EXAMPLE  PROBLEM 

A  small  three  period  example  is  considered  in  this  section. 

The  basic  network,  consisting  of  four  nodes  and  ten  arcs,  is  illustrated 
in  Figure  3.  Three  of  the  nodes  have  loop  arcs  associated  with  them 
to  capture  their  storage  capability.  The  arc  data  for  the  basic  net¬ 
work  is  summarized  in  Table  1,  and  the  net  demands  for  each  node  are 
presented  in  Table  2.  Thirteen  units  are  available  in  the  distribution 
system  during  period  one.  Eight  units  are  available  during  period  two, 
and  four  units  are  available  during  period  three. 


FIGURE  3 
BASIC  NETWORK 
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The  matrix  formulation  of  the  capacitated  transshipment 
problem  is  given  in  Figure  4,  and  the  corresponding  master  network 


is  illustrated  in  Figure  5.  This  figure  highlights  some  of  the 
important  topological  characteristics  of  the  multi-period  trans¬ 
shipment  problem.  The  nodes  of  the  master  network  are  arranged  in 
m  parallel  rows  and  T  +  1  parallel  columns.  Row  v  of  the  master 

network  corresponds  to  node  v  of  the  basic  network.  Column  t  of 

th 

the  master  network  corresponds  to  the  t  point  in  time,  or  the  end 
of  period  t.  Figure  5  illustrates  that  the  master  network  is  an 
acyclic  graph.  All  master  arcs  are  directed  from  a  node  in  column 
t  -  1  to  a  node  in  column  t.  The  movement  of  units  in  the  basic 
network  during  period  t  corresponds  to  the  flows  on  the  arcs  directed 
from  column  t  -  1  nodes  to  column  t  nodes.  The  pattern  of  arcs  be¬ 
tween  columns  t  -  1  and  t  is  repeated  for  each  period.  The  marginal 
costs  and  upper  bounds  on  these  arcs  are  also  identical  from  period 
to  period. 

Figure  6  shows  a  feasible  solution  to  the  master  network 
problem,  and  Figure  7  shows  the  corresponding  solution  for  the  basic 
network  for  each  period.  This  solution  has  an  objective  function 


value  of  414. 


FIGURE  5 
MASTER  NETWORK 


6.  PNET-MP:  A  MULTI-PERIOD  TRANSSHIPMENT  CODE 


A  specialized  version  of  the  efficient  primal  simplex  network 
code  PNET  [1]  is  presented  in  this  section.  This  large-scale  implementa¬ 
tion  uses  the  same  principal  list  structures  as  PNET  to  store  the 
tree  representation  of  the  basis  matrix.  Specifically,  the  M  node 
and  M  -  1  arc  triangular  basis  matrix  is  represented  by  six  node 
functions. 

The  underlying  structure  of  the  basis  tree  is  described  by 
the  predecessor  node  and  predecessor  arc  functions.  The  thread  and 
depth  functions  are  used  to  improve  the  pivoting  capabilities  of  the 
code.  The  primal  basic  solution  and  the  associated  dual  solution  are 
also  maintained  as  node  functions.  Each  of  these  node  functions  is 
implemented  in  PNET-MP  as  a  node  length  array. 

A  basis  tree  corresponding  to  the  feasible  master  network 
solution  from  Figure  6  is  shown  in  Figure  8.  Table  3  provides  the 
six  node  functions  that  are  used  to  describe  the  basis  tree  in  Figure 
8. 

The  principal  difference  between  the  general  purpose  code, 

PNET,  and  the  specialized  multi-period  code,  PNET-MP,  involves  the 
method  of  representing  the  original  problem  data.  PNET  uses  three  N 
length  lists  to  store  the  arc  data.  That  is,  the  destination  node, 
marginal  cost,  and  upper  bound  of  every  master  arc  is  explicitly 
stored  in  the  arc  lists.  This  means  that  each  basic  network  arc  is 
duplicated  T  times  in  the  arc  lists.  PNET-MP,  on  the  other  hand, 
uses  two  n  length  and  one  m  length  lists  to  store  the  arc  data  of 
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TABLE  3 

NODE  FUNCTIONS 


I 

PN(I) 

PA(I) 

THREAD (I) 

DEPTH (I) 

NET CAP (I) 

POT (I) 

1 

0 

0 

6 

0 

0 

0 

2 

5 

4 

7 

2 

oo  -1 

-7 

3 

7 

7 

8 

2 

4  -1 

28 

4 

6 

9 

5 

2 

00  -0 

-28 

5 

1 

1 

2 

1 

3 

0 

6 

1 

2 

9 

1 

0 

10 

7 

1 

3 

11 

1 

6 

30 

8 

3 

8 

1 

3 

2 

40 

9 

6 

14 

13 

2 

1 

17 

10 

7 

16 

3 

2 

5 

55 

11 

7 

17 

16 

2 

2 

32 

12 

14 

29 

15 

4 

oo_l 

19 

13 

9 

21 

4 

3 

1 

17 

14 

11 

26 

12 

3 

1 

57 

15 

11 

27 

10 

3 

1 

34 

16 

11 

28 

14 

3 

0 
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the  basic  network.  The  N  master  network  arcs  are  implicitly  stored 
by  only  explicitly  storing  the  n  basic  network  arcs. 

Since  only  one  copy  of  each  basic  network  arc  is  physi¬ 
cally  stored,  PNET-MP  has  the  capability  to  solve  extremely  large 
multi-period  network  problems.  For  instance,  PNET  requires  over 
two  and  a  half  times  more  core  storage  than  PNET-MP  to  solve  a 
five  period  problem  with  a  100  node  and  500  arc  basic  network. 

For  a  ten  period  problem  with  a  100  node  and  1000  arc  basic  net¬ 
work,  PNET  requires  over  six  times  more  core  storage  than  PNET-MP. 
As  either  the  number  of  periods  or  the  arc  to  node  ratio  increases, 
PNET-MP  becomes  more  attractive  than  PNET  in  terms  of  total  storage 
requirements. 

As  pointed  out  in  Section  4,  each  column  of  the  M  x  N 
constraint  matrix  corresponds  to  a  master  arc  and  each  row  corres¬ 
ponds  to  a  master  node.  Master  arc  K,  from  master  node  I  to 

tv 

master  node  is  related  to  arc  k  of  the  basic  network  for  period 
t  in  the  following  manner: 

K  *  k  +  (t  -  l)n 

'k  '  \  +  (t  -  1)0 

JK  ■  Jk  +  “ 

PNET-MP  uses  this  fundamental  relationship  between  the  master  and 
basic  arcs  to  generate  the  complete  master  network  even  though 
only  the  basic  network  is  explicitly  stored. 

The  arcs  of  the  basic  network  are  stored  in  two  n 
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length  and  one  m  length  lists.  The  destination  node  of  basic 
arc  k  is  stored  in  TO(k).  The  arcs  are  stored  according  to 
common  origin  node.  The  destination  node  of  the  first  arc  out 
of  each  origin  node  is  negated  in  the  TO  list.  This  is  used  to 
signal  the  start  of  the  arc  data  for  the  next  node.  The  mar¬ 
ginal  cost  of  basic  arc  k  is  stored  as  COST(k).  The  upper  bound 
on  the  loop  arc  associated  with  basic  node  v  is  stored  as  CAP(v). 
This  node  length  list  is  used  to  represent  the  finite  arc  capa¬ 


cities.  That  is,  basic  arc  k  has  an  infinite  upper  bound  if 

i“k  f  and  has  a  finite  upper  bound  of  CAP(i^)  if  i^  *  Table 

4  summarizes  the  arc  data  for  the  four  node  and  ten  arc  basic  net¬ 


work  given  in  Section  5.  Note  that  a  trailer  record  (position  11) 


is  used  to  signal  the  end  of  the  arc  data. 


TABLE  4 
ARC  FUNCTIONS 


This  particular  data  structure  is  well-suited  for  sequential 
processing  of  the  arc  data.  Extensive  computational  testing  has  shown 
that  sequential  processing  is  highly  effective  for  primal  simplex  net¬ 
work  optimization  codes  [10,  11,  16,  17]. 

There  are  four  fundamental  operations  that  must  be  performed 
by  a  primal  simplex  network  code.  They  are: 

1.  Determination  of  an  initial  solution, 

2.  Selection  of  an  entering  arc, 

3.  Selection  of  a  leaving  arc,  and 

4.  Pivot,  or  change  of  basis. 

The  precise  implementation  in  PNET-MP  of  each  of  these  operations  will 
be  considered  next. 

An  initial  basis  for  the  problem  can  be  obtained  by  adding 

an  artificial  variable  for  each  node  constraint  of  the  master  network. 

th 

The  coefficient  of  the  artificial  variable  in  the  I  constraint  is 
determined  by  the  sign  of  the  net  demand  for  units  at  node  I.  If  node 
I  has  a  non-negative  net  demand,  then  the  artificial  variable  has  a 
+1  coefficient  in  constraint  I.  Otherwise,  it  has  a  -1  coefficient. 
Each  artificial  variable  is  assigned  an  infinite  marginal  cost.  The 
expanded  master  network  problem  is  given  by: 
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Minimis.  CX1  +  CX2  *  CX3. . .  +  CX1*1  +  CXT  *•'  X  H°  *■  X  X1  *  M  V2 


+  n  *T“l  f  iir 


iubjtce  co: 


ix1  >  rr 


♦  W1 


4-  D*y  -  b  ‘ 


0  S  XC  £  8  for  6  •  l.  2,  .. 
0  £  YC  cor  c  »  0,  1,  .. 


where  M  is  an  m  dimensional  infinite  (Big  M)  vector,  and  D  is  a 

th 

diagonal  matrix  whose  v  diagonal  coefficient  is  given  by 

„  1  if  bC  >  0 

*  v 

w  -1  if  bc  <  0 

V 

This  expanded  problem  has  full. row  rank  since  each  D* 
has  rank  m.  An  equivalent  expanded  pure  network  problem  can  be 
obtained  by  creating  a  redundant  constraint  equal  to  the  negative 
of  the  sum  of  the  M  network  constraints.  This  new  constraint  has 


the  form 


d°Y°  +  d1^  +  d2Y2*..  +  d1"^1”1  +  dTYT  -  0 


where  dc  is  an  m  dimensional  vector  whose  vC^  coefficient  is  the 
negative  of  the  vC^  diagonal  coefficient  of  DC.  This  redundant 
constraint  can  be  thought  of  as  an  artificial  node  in  the  ex¬ 
panded  master  network.  The  M  artificial  variables  correspond  to 
arcs  directed  into  and  out  of  the  artificial  node.  Since  Xfc  ■  0 


and  YC  *  Dtbt  is  a  feasible  solution  to  the  expanded  network  problem, 
the  M  artificial  arcs  can  be  used  as  an  initial  basis  tree.  For 
convenience,  the  artificial  node,  M  +  1,  is  used  as  the  root  of  the 
basis  tree.  Figure  9  illustrates  this  starting  basis  tree  for  the 
small  example  in  Section  5. 

PNET-MP  uses  this  artificial  arc  basis  tree  as  the  ini¬ 
tial  feasible  solution.  The  six  node  functions  are  extremely  easy 
to  initialize  for  this  starting  basis. 

A  modification  of  the  outward  node  most  negative  pivot 
rule  [11]  is  used  by  PNET-MP  to  select  the  arcs  to  enter  the  basis. 
The  modification  is  motivated  by  the  fact  that  each  arc  of  the 
basic  network  generates  T  master  arcs.  Instead  of  selecting  the 
most  pivot  eligible  arc  out  of  a  single  master  node,  PNET-MP 
selects  the  best  arc  out  of  all  T  related  master  nodes.  That  is, 
the  arcs  originating  at  master  nodes  i,  i  +  m,  i  +  2m,  ...,  and 
i  +  (T  -  l)m  are  examined  together. 

PNET-MP  uses  the  standard  upper  bounded  simplex  minimum 
ratio  test  to  determine  the  leaving  arc.  This  is  carried  out  by 
determining  the  maximum  allowable  flow  change  on  the  unique  cycle, 
or  basis  equivalent  path,  created  by  adding  the  entering  arc  to 
the  current  basis  tree.  If  more  than  one  arc  reaches  a  bound, 
then  the  first  such  arc  identified  is  selected  to  leave  the  basis. 

A  slightly  more  complicated  selection  rule  has  been  developed  [  2 , 

4  ,  7  ,  8  ]  that  guarantees  the  finiteness  of  the  network  simplex 
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algorithm,  but  it  was  not  implemented  in  order  to  simplify  com¬ 
parisons  with  PNET,  which  uses  the  standard  rule. 

PNET  uses  a  simple  flag  system  to  identify  the  arcs  that 
are  non-basic  at  upper  bound.  This  system  involves  the  negation  of 
the  arc  length  capacity  array.  PNET-MP  cannot  use  this  method 
since  only  the  upper  bounds  of  the  loop  arcs  of  the  basic  network 
are  explicitly  stored.  Instead,  PNET-MP  uses  a  simple  bit  map 
technique  to  indicate  the  non-basic  status  of  the  master  network 
arcs.  Due  to  the  special  structure  of  this  multi-period  problem, 
the  bit  map  can  generally  be  implemented  quite  easily  with  only  a 
single  m  length  array. 

PNET-MP  performs  the  actual  basis  exchange  operation,  or 
pivot,  in  a  manner  identical  to  that  used  in  PNET.  The  necessary 
updates  of  the  various  node  functions  are  described  in  detail  in 
[  5,  12,  14]. 

In  order  to  properly  assess  the  performance  of  PNET-MP, 
twelve  large  multi-period  test  problems  were  generated.  Table  5 
shows  the  dimensions  of  both  the  basic  and  the  master  networks  for 
each  of  the  test  problems.  The  table  also  presents  the  computation¬ 
al  results  on  these  problems  for  the  specialized  code,  PNET-MP,  as 
well  as  the  general  purpose  pure  network  flow  code,  PNET. 

In  terms  of  core  storage  requirements,  PNET-MP  is  clearly 
superior  to  PNET  for  this  class  of  problems.  For  Instance,  on  the 
largest  problem  in  the  test  set,  PNET  requires  more  than  four  times 
the  space  of  PNET-MP.  In  fact,  due  to  excessive  core  requirements. 


In  60,000  words 


PNET  was  unable  to  solve  three  of  the  problems  on  The  University  of 
Texas1  CDC  6600. 

Solution  statistics  reported  in  Table  5  include  the  number 
of  pivots,  optimization  time,  solution  time,  and  mean  time  per  pivot. 
All  reported  times  are  in  c.p.u.  seconds  on  the  CDC  6600  using  the 
MNF  FORTRAN  compiler.  Times  do  not  include  any  input  or  output 
operations.  Solution  time  is  equal  to  the  sum  of  the  start  time 
(initialization  and  construction  of  the  first  basis  tree)  and  the 
optimization  time.  The  mean  time  per  pivot  is  simply  the  ratio  of 
the  optimization  time  to  the  number  of  pivots. 

At  first  it  was  surprising  that  PNET-MP  actually  solved 
the  problems  faster  than  PNET  since  PNET-MP  was  designed  primarily 
to  reduce  core  storage,  not  to  execute  faster.  Upon  closer  ex¬ 
amination  it  was  determined  that  the  success  of  PNET-MP  revolves 
around  the  pivot  strategy  that  it  employs.  Specifically,  PNET-MP 
processes  all  of  the  master  network  arcs  associated  with  a  single 
basic  network  arc  together,  whereas  PNET  processes  these  master 
network  arcs  separately.  As  indicated  by  the  number  of  pivots, 
the  pivot  strategy  used  by  PNET-MP  is  superior,  for  this  class  of 
problems,  to  that  used  by  PNET.  As  expected,  PNET-MP  requires 
more  time  per  pivot  than  PNET  since  it  must  perform  more  work  in 
order  to  construct  the  master  network  from  its  single  copy  of  the 


basic  network. 
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